
% laod boot results:
load(BootResults);
theta_est  = theta_boot; % below only reference to theta_est.
%pproc      = pproc_boot; % replace pproc with bootstrapped processes.
p_unob_het = p_uh_old; 
% load corresponding options:
load(EstimatesResults, 'opt', 'pproc');

%% Compute long-run of current process:
% Prepare data, states, pproc etc based on options from result file 
% (typically result file will not include this data structure)
optCount = opt; optCount.SelSubSample = 0; 
optCount.TrunkAge = 0; % Option to trucate age (added later, so may not be present in some results files)
[FID_Grid,Action,Age,AgeState,States,~,DS] = prepdataM10(optCount);

% Get Baseline value function
[V,Pa,F] = ValueFunction(States,pproc,theta_est,opt);
% Create boxes for sugarcane proximity:
BoxSize = 2000*States.SRDMidPoints(1:end-1,1)';
[Boxes] = BoxesCreate([DS.CoordX,DS.CoordY], BoxSize);

[UHprob,~] = UpdateUHProb2(p_unob_het,States,AgeState,Action,Pa);

IniAgeStates = double(AgeState(:,end)); IniExoTimeStates = double(States.ExoTimeStates(:,end));

Ns = 80;

CF.FID_Grid = FID_Grid;

%disp(' Running simulations for baseline case: ')

for s=1:Ns
    [ OutputPathBase(:,s),  YieldPathBase(:,s), AcreagePathBase(:,s), LandUseBase(:,s), OutputDetailBase(:,s)]  = SimulationPathsLong(IniAgeStates, IniExoTimeStates, Pa, States,Boxes, pproc, theta_est, UHprob, opt );
    %disp(' ')
    %ProgressBar(s,Ns)
end

CF.OutputPathBase = OutputPathBase; CF.YieldPathBase = YieldPathBase; 
CF.AcreagePathBase = AcreagePathBase; CF.LandUseBase = LandUseBase; 
CF.OutputDetailBase = OutputDetailBase;

CF.LandUseBase = mean(CF.LandUseBase,2);
CF.OutputDetailBase = mean(CF.OutputDetailBase,2);

%fprintf('\n');

save(SaveCFto, 'CF', '-v7.3');


%% Price change counterfactual
% Get States and pproc for changed price process:
PriceChange.on = 1;
PriceChange.factor = 1.05;
[StatesChange,pprocChange] = GenStateVariables(DS,optCount,double(AgeState),PriceChange);
% Get price change value function
[VChange,PaChange,FChange] = ValueFunction(StatesChange,pprocChange,theta_est,optCount);

for s=1:Ns
    [ OutputPathCF05(:,s),  YieldPathCF05(:,s), AcreagePathCF05(:,s), LandUseCF05(:,s), OutputDetailCF05(:,s)]  = SimulationPathsLong(IniAgeStates, IniExoTimeStates, PaChange, StatesChange, Boxes, pprocChange, theta_est, UHprob, opt );
    %disp(' ')
    %ProgressBar(s,Ns)
end

CF.OutputPathCF05 = OutputPathCF05; CF.YieldPathCF05 = YieldPathCF05; 
CF.AcreagePathCF05 = AcreagePathCF05; CF.LandUseCF05 = LandUseCF05; 
CF.OutputDetailCF05 = OutputDetailCF05;

CF.LandUseCF05 = mean(CF.LandUseCF05,2);
CF.OutputDetailCF05 = mean(CF.OutputDetailCF05,2);

%fprintf('\n');

save(SaveCFto, 'CF', '-v7.3');
